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Abstract 

We study the bond percolation problem in random graphs of N weighted vertices, where 
each vertex i has a prescribed weight Pi and an edge can connect vertices i and j with 
rate PjPj. The problem is solved by the q — > 1 limit of the g-state Potts model with inho- 
mogeneous interactions for all pairs of spins. We apply this approach to the static model 
having Pi oc i~f(0 < /u < 1) so that the resulting graph is scale-free with the degree exponent 
X = 1 + l//u. The number of loops as well as the giant cluster size and the mean cluster size 
are obtained in the thermodynamic limit as a function of the edge density, and their associ- 
ated critical exponents are also obtained. Finite-size scaling behaviors are derived using the 
largest cluster size in the critical regime, which is calculated from the cluster size distribu- 
tion, and checked against numerical simulation results. We find that the process of forming 
the giant cluster is qualitatively different between the cases of X > 3 and 2 < X < 3. While 
for the former, the giant cluster forms abruptly at the percolation transition, for the latter, 
however, the formation of the giant cluster is gradual and the mean cluster size for finite N 
shows double peaks. 



1 Introduction 



In the last few years graph theoretic approach has been of great value to charac- 
terize complex systems found in social, informational and biological areas. Here, 
a complex system is represented as a graph or network whose vertices and edges 
stand for its constituents and interactions. A simplest model for such is the random 
graph model proposed by Erdos and Renyi (ER) [1]. In the ER model, N num- 
ber of vertices are present from the beginning and edges are added one by one 
in the system, connecting pairs of vertices selected randomly. Due to the random- 
ness, the distribution of the number of edges incident on each vertex, called the 
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degree distribution, is Poissonian. However, many real-world networks such as the 
World-wide web, the Internet, the coauthorship, the protein interaction networks 
and so on display power-law behaviors in the degree distribution. Such networks 
are called scale-free (SF) networks [2]. Thanks to recent extensive studies of SF 
networks, various properties of SF network structures have been uncovered [3,4,5]. 

There have been a few attempts to describe scale-free networks in the framework of 
equilibrium statistical physics, even though the number of vertices grows with time 
in many real- world networks [6,7,8,9,10]. In this approach, various mathematical 
tools developed in equilibrium statistical physics may be used to understand net- 
work structures. To proceed, one needs to define equilibrium network ensembles 
with appropriate weights, where one graph corresponds to one state of the ensem- 
ble. In a canonical ensemble, the number of edges L is fixed: Given a degree dis- 
tribution, Pd(k), the mean degree (k) = Y*kpd(k) is obtained. Then the number of 
edges obtained through the relation, L = (k)N/2, can be fixed. A degree sequence 
specifies the number of vertices with degree k as Pd(k)N [11,12]. 

A grandcanonical ensemble can be also defined, where the number of edges is 
also a fluctuating variable while keeping the SF nature of the degree distributions. 
The grandcanonical ensemble for SF random graphs is realized in the static model 
introduced by Goh et al. [13] or in its generalized version investigated in Refs. [14]. 
The name 'static' originates from the fact that the number of vertices is fixed from 
the beginning. Here each vertex i has a prescribed weight Pi summed to 1 and an 
edge can connect vertices i and j with rate PjPj [15,16,17]. A chemical potential- 
like parameter K that can be regarded as "time" in the process of attaching edges 
controls the mean number of edges so that (L) increases with increasing K. 

As the parameter K increases, a giant cluster, or giant component, forms in the sys- 
tem. Here the giant cluster means the largest cluster of connected vertices whose 
size is O (N) . Often such a giant cluster appears at the percolation transition point. 
In equilibrium statistical physics, the percolation problem can be studied through 
a spin model, the g-state Potts model in the q — > 1-limit [18]. Using the relation, 
in this paper, we study the evolution of SF random graphs from the perspective of 
equilibrium statistical physics. To be specific, we construct the g-state Potts model, 
where the interaction strength between each pair of vertices is inhomogeneous on 
the complete graph. In this formulation, since the interaction strength K is tunable, 
the mean number of edges (L) varies. Thus, the grandcanonical ensemble is taken 
in the network representation. However, since the number of spins (vertices) is 
fixed, the formulation corresponds to a canonical ensemble in the spin-model rep- 
resentation. Note that our model is different from the one studied by Dorogovtsev 
et al, [19] where the Potts model is defined as a given fixed network so that each 
edge represents homogeneous interactions. 

The formulation of the spin model facilitates explicit derivation of various prop- 
erties of the SF network. Thus we derive the formula for the giant cluster size, 
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the mean cluster size, and in particular, the number of loops and clusters. These 
quantities are explicitly evaluated analytically for the static model with 
(0 < fx < 1) in the thermodynamic limit as a function of the edge density, and 
their critical properties are also studied. The degree exponent X, is related to by 
X, = 1 + 1///. Moreover, their finite-size scaling behaviors are obtained using the 
finite largest cluster size for finite N that in turn is evaluated from the cluster size 
distribution. From these, we are able to elucidate the process of formation of the 
giant cluster. While for the case X, > 3, the giant cluster forms abruptly at the per- 
colation transition point K c , for the case 2 < X, < 3 where most real world networks 
belong to, however, the formation of the giant cluster is gradual and the mean clus- 
ter size for finite N show double peaks. 

In fact, the percolation problem of SF networks has been studied, but in a different 
way, that is, by removing vertices one by one as well as their attached edges from an 
existing SF network [20,21,22]. The percolation transition was understood by using 
the branching process approach, which is supposed to be valid near the percolation 
transition point, where the network is sparse. In this paper, we provide the criterion 
for the validity of the branching process approach for a general degree distribution, 
and show that the branching process and the Potts model approaches are equivalent 
for the static model. Finally, note that while the branching process approach cannot 
count the number of loops, the Potts model formalism we use here enables us to 
count it. 

This paper is organized as follows. We introduce in Sec. 2 an ensemble of ran- 
dom graphs where each vertex is weighted, and present in Sec. 3 the Potts model 
formulation to derive graph theoretical quantities from its free energy. In Sec. 4, 
the connection between the Potts model formulation and the branching process ap- 
proach is discussed. The general results of Sec. 3 are applied to the static model 
in Sec. 5 to obtain explicitly the giant cluster size, the mean cluster size and the 
mean number of loops and clusters as a function of K. The cluster size distribution 
and the largest cluster size in finite size systems are obtained in Sec. 6. The finite- 
size scaling is presented and compared with numerical simulation results in Sec. 7. 
Finally Sec. 8 contains summary and discussion. 



2 Random graphs with weighted vertices 

Suppose that the number of vertices /V is fixed (static) and each vertex i = 1 , . . . ,N is 
given a probability P, summed to 1 . The ER model of random graphs corresponds 
to assigning Pi = l/N for all i. To construct a SF graph, we use Pi ~ i~ M with 
< fi < 1 . However, for the time being, Pi is arbitrary as long as Pi -C 1 for all i. 

In each unit time duration, two vertices i and j are selected with probabilities P, and 
Pj. If i = j or an edge connecting i and j already exists, do nothing; otherwise, an 



3 



edge is added between the vertices i and j. This process is repeated for times. 
Then the probability that a given pair of vertices i and j (i ^ j) is not connected by 
an edge is given by (1 - 2PiPj) NK ~ e - 2NKp i p j, while that it does is 1 - e - 2NKp i p j, 
Here we used the condition P,- <C 1 . The factor 2 comes from the equivalence of 
(/y) and (ji). We use the "interaction" parameter K for later convenience which 
controls the edge density (L)/N. Since each edge by is produced independently, 
this process generates a graph G with probability 



P(G)= (\-e- 2NKPiP J) ]^[ e~ 2NKP ^ 



b u €G 

-NK{\-M 2 ) j-j (JNKPlPj _ ^ (1) 



e 



b u €G 



where we used the notation M n = J^lj P". By a graph G, we mean a configuration 
of undirected edges connecting a subset of N(N — l)/2 pairs of labelled vertices 
i= 1,2,..., AT. 

We then evaluate the ensemble average of any graph theoretical quantity A by 

(A)=£P(G)A(G). (2) 

G 

One example is the degree /c ( - of a vertex z, the number of edges incident on i. To do 
this, the generating function of /c ( -, g,-((fl) = (co^'), is first calculated as 

g ( -(to) = [e- 1NKP ' P > + to( 1 - • (3) 

j(/0 



From this, one has 

= £(l- e - 2 ^'), (4) 



and the average degree (k) is 



Also, 
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<*?>=(«>£) „(») 



(0=1 



j;(i- e -^) 

= <**>+&> 2 - 



+ £ e -2NKPiP j ( l _ e -2NKP i Pjj 



-IN KPjPj \ 2 



(6) 



We remark that Eq. (3) is rewritten as 

g f (eo) ~ e -( 1 -«»)(*'-> 



(7) 



with (ki) in Eq. (4) when |( 1 — co) ( 1 — e - 2NKPiP i) | < 1 for all /) . It implies that 
the probability that is equal to k, pd,i(k) = (§kj,k), is given by 



Pd,i(k) 



1 d fc 



(0=0 



(8) 



for k 3> 1. Other quantities are discussed later on. 



3 Potts model 



3. 1 Potts model and random graph 



It is well known that the g-state Potts model provides a useful connection between 
the geometric bond percolation problem and the thermal systems through the Kaste- 
leyn construction [18]. The q — > 1 limit of the Potts model corresponds to the bond 
percolation problem. The same approach can be used for the random graph prob- 
lem. From the viewpoint of the thermal spin system, this is basically the infinite 
range model since all pairs of spins interact with each other albeit with inhomoge- 
neous interaction strength. 

Consider the g-state Potts Hamiltonian given by 

N 

-H = WK^PiPjHohOj) +/*oI>8(g ! -, 1) - 1], (9) 

i>j i=l 

where K is the interaction, ho is a symmetry-breaking field, 8(x,y) the Kronecker 
delta function, and o ; the Potts spins taking integer values 1,2, ...,q = r+1. We 
use the notation r = q — 1. The partition function Z N (q 7 ho) can be written as 
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Z N (qM) =Tre H = Trf[ [l + {e 2NKP ^ - l)5(o ( -,o,-)] 

xJ-JgAoWo,-,!)-!)^ (1Q) 



where Tr denotes the sum over q N spin states. Expanding the first product and 
taking the Tr operation, one has 

z N {qM) =e n (^^-i)n(^+ re " /!o ) nG( " ) ' en) 

G 2> (7 €G s>l 

where nc(s) is the number of s-clusters, a cluster with s vertices in a given graph 
G. Comparing this with Eq. (1), one immediately notices that 

Z N (q,h ) = e NK(l - M ^ £P(G) Yl(e srh ° + re- sh °) n °W. (12) 

G s>l 

In particular, Z N (q,0) = e NK ( l ~ Ml ) (q c ), where C = L^g^) is the total number of 
clusters in graph G. Thus Z/v(g,0) is the generating function of C. 

The magnetization of the Potts model is defined as 

m(q,ho) = ^^-\nZ N (q,h ). (13) 

It can be written as 

e irh Q+re -sh 



, ^clusters 

m(q,h ) = 



(§)T\ c iusUe srhQ + re- shQ ) 

(14) 



(nclusters(^° + ^-^)) 



If we introduce the cluster size distribution P(s) = n(s)(s/N) with n(s) = (nc(s)) 
and the generating function <P (z) = 'Ls>iP( s )z s , the magnetization is, when q=l, 

m{l,h ) = U-e- sh A = \-?(e- h °). (15) 

s>l V J 

The generating function fp(z) will be used in Sec. 6 to investigate the asymptotic 
behavior of the cluster size distribution. 

When ho = 0, the magnetization vanishes for finite N. However, when we take the 
limit ho — > after the thermodynamic limit N — > °°, the contribution from the largest 
cluster whose size is S can survive to give 



m(l,ho — > 0) = 




if S//V is finite. Let us define a giant cluster by a cluster whose size is (N). Then 
m(l,/?o — > 0) is the ratio of the giant cluster size to A 7 , if it exists, and the system 



6 



q=2 q=3 q-4 



t 



Fig. 1. Vector representations of g-state Potts spins with q = 2, 3, and 4. 

is considered as being in the percolating phase if m(l,/zo — > 0) is non-zero. For 
simplicity, we will call m(l,ho — > 0) the giant cluster size and denote it by m. 

The susceptibility defined as %{q,ho) = (l/q)(d/dho)m(q,ho) on the other hand is 
related to the mean cluster size: 

%(l,h - 0) = lim lim J>( S > e -^ = £ sP( s ), (17) 



where /lo(S') — > °° is used with (S) the ensemble average of the largest cluster size, 
which we call simply the largest cluster size. We will denote %( 1 , ho — > 0) by s. Note 
that our definition of s is normalized with respect to the total number of vertices 
instead of the number of vertices belonging to finite clusters. 

The number of loops N\ oov is related to the total number of clusters C through 

N loop = L-N + C. (18) 

Since Z N (q,0) ~ (g c ),one can notice that the number of loops per vertex (N\ oop )/N 
is given as 

^ = ^-l + i |-[lnZ^,0)L !• (19) 
N N Ndq 1 W Jlq - 1 

We will denote (Moop) /N by £ and call it the number of loops for simplicity. 



3.2 Partition function 



A convenient way to evaluate the partition function is to resort to the vector-spin 
representation where one associates an r-dimensional vector 5(o;) of unit length to 
each spin value a,-, where 5(1) = ( 1 , 0, . . . , 0) and 5(o,-) with a,- = 2, 3, . . . , q point 
to the remaining r corners of the r-dimensional tetrahedron (See Fig. 1). Then one 

— * 

can represent the Kronecker delta function as a dot product between S's, 

d(c i ,Cj) = -(rS(c i )-S(c j ) + l). (20) 
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Using this, the interaction term in Eq. (9) can be written as 

INKj^PiPfiia, Oj) = NK Q -M^j + ^ (^PiS(Oi) j . (21) 

The perfect square is then linearized through the identity / dye~ ay2+by = ^n/ae b2 ^ 4a ^ . 
Thus we have 



Z N (q 1 h ) =Tr 



(22) 

where the integration is over the r-dimensional space and ho = (ho, 0, . . . , 0). Now 
the Tr operation can be performed for each spin independently. Defining 



1(h) = - £ (23) 

<J=1 

one then has 

Z N (q,ho) = {^) L \ N ^ Ml) J dye-^+ZlMi^m. (24 ) 



Provided (\/N)^ =l \nl(h + NPiy) has a well-defined limit as N —> °°, one can 
apply the saddle point method to Eq. (24), where the integral is replaced by the 
value of the integrand at its maximum. The maximum is for y which is a solution 
of the saddle-point equation 

^y=^t V M(ho + NPiy). (25) 

When ho = 0, y = is always a solution of Eq. (25), but the spontaneous symmetry- 
breaking solutions (y ^ 0) with the Potts symmetry may appear for large K. When 
the symmetry-breaking field is applied along the 1 -direction, ho = (ho, 0, . . . , 0), the 
nontrivial physically relevant solution of Eq. (25) is expected in the sub-manifold 
y = (y, 0, . . . , 0) . The limit ho — > + then selects one of the q equivalent spontaneous 
symmetry -breaking solutions. With this in mind, we may restrict our attention to 
the one-dimensional sub-manifold of y in Eqs. (24) and (25). As a result, we then 
have, as N — > °°, 

^lnZ N (q,h ) = Xnq + K^-M^j -rF(y,h ), (26) 

with 

F{yM) = ^y 2 -^t^(ho + NP iy ,q), (27) 
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where 

«*,) = !£«*<■> = £?±£l, (28) 
and y is the solution of the one-dimensional saddle-point equation, 

&=\% p >k XaUh>+NP ' y > q) - (29) 

Here, the (/-dependence in F(y, ho) is not shown explicitly. Since m = ( 1/ (rN))(d/dho) InZ, 
we see that 

m{qM) = -^bM = ^j^^ln^ho + NP^q). (30) 

At this point, it is useful to take the r — > limit in Eqs. (27), (29), and (30), which 
yields, with hi = ho + NPty, 

F(y,h ) = ^J 2 ~ ^ J> _ * ~ 1 (3D 



where y is the solution of 



AIT Nf- =x 



N 



^ = Y,m-e- hi ). (32) 
(=l 

The magnetization and the susceptibility reduce to 

m(l,*o) = -l£(l-«-**), (33) 

and 



X(l./.o)4|^(l + ^). 

= if e ->, + (a,*-*) 2 (34) 



respectively, where it is used that 

Y N . ft*"** 

(35) 



Thus the giant cluster size and the mean cluster size are obtained from Eqs. (33) 
and (34), respectively, with ho — > 0. 
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Also, the number of clusters per vertex is 



{ ^- = l-K-F(y,0)\ q=h (36) 



while that of loops is 

(Moop) _ (L) 
N N 

where y is given by Eq. (32) with ho = 0. 



K-F(y,0)\ q=h (37) 



When ho — > 0, a non-trivial solution of Eq. (32) begins to appear when (2K) 1 < 
N^f =l P 2 , which gives the following characteristic value K c 

K c = 1 — ~ • (38) 

2A^f =1 /f 

When decays slower than z' _1 and <C 1 for all i, (k) = 2K and (k 2 ) = 
(l/N)^f =l (kf) = (k) +N- l ^ l (k i } 2 = 2K + 4NK 2 ZtiPh which will be shown 
below. Then the condition K = K c is equivalent to the well-known condition (k 2 ) / (k) 
2 [11]. Whether the percolation transition occurs at K c or not will be investigated 
for specific P^s of the static model. 



4 Branching process approach 



The cluster size distribution P(s) can be obtained from Eqs. (32) and (33) using 
m(l,ho) = l—^p(e~ h °) andi^s) = (l/s\)d s P (z) /dz s \ z =o- However, the cluster size 
distribution can also be obtained through the generating function approach or equiv- 
alently the branching process approach. Here, the presence of loops infinite clusters 
is neglected and each cluster in a given graph is considered as a tree generated by 
successive branchings from an arbitrary vertex [23,24]. Consider the probability 
that a randomly chosen vertex belongs to a ^-cluster, which is just P(s). Then P(s) 
can be written recursively as 

P(s) = 8 MW (0) + £ p d (k) n£*(*,) S L .,„-i, (39) 

k>\ i=l Si 

where Pd{k) is the degree distribution and R(s) is the probability that a randomly- 
chosen edge has a ^-cluster at its one end, and thus equal to the number of edges 
followed by ^-clusters divided by 2L. R(s) is obtained self-consistently as 

R(s) = 5.,!r d (0) + £ r d (k) fl^R( Sl ) 5 L .,,,_!, (40) 

k>\ i=\ si 
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where r^(k) is the probability that the vertex at either end of a randomly-chosen 
edge has k + 1 edges and thus is equal to (k+l)pd(k+l)/{k). With the generating 
functions <p(z) = I£ =l P(s)z s and ^(z) = l^ =l R(s)z s , Eqs. (39) and (40) can be 
written in more compact forms as 

v(z)=zg(lL(z)) (41) 

and 

2L(z) = zfMz)), (42) 

where g(a>) = I£ = oP d (k)rf and /(to) = g'(co)/(&> = g'(co)/g'(l). We mention 
that Eqs. (41) and (42) with z = 1 are equivalent to those derived by Molloy and 
Reed [1 1] for a given degree sequence. 

Eqs. (41) and (42) are the standard results. For the grandcanonical ensemble we are 
using, the generating functions g(cfl) and /(co) are represented in terms of gi((o), 
the generating function of Pd,i(k) in Eq. (3), as 



1 N 

iV (=1 

/((D) = J-f^ (43) 
J{ ' N(k}L d& ^ 

with (k) = g'(l). In particular, if 

^■((D) = e^ 1 -^™ (44) 

which holds, for example, when 1 — e NKP ' p j <C 1 for all i ^ j, then Eqs. (41) and 
(42) of the branching process approach are exactly equal to Eqs. (33) and (32) of 
the Potts model formulation, identifying z, 2 (z), and Hi (z) with e~ h °, 1 — m(l, ho), 
and l—y/ {IK), respectively. 



5 Percolation of the Static model: Thermodynamic limit 



So far, our discussion applies to arbitrary In this section, we specialize to the 
case of the static model, 

Pi = rT^ (0 </*<!)• (45) 
wW 

Here C,n(x) = Y4L1 i~ x > an d in the limit N — > °°, it converges to the Riemann zeta 
function C,(x) when x > 1 and diverges as A^ 1_x /(1 — x) when < < 1. In the 
marginal case with x = 1, it is calculated as ln/V + yM with y^ = 0.5772..., the 
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Fig. 2. Giant cluster size m, mean cluster size s, and number of loops I = {Ni oov )/N vs. 
K for n = 5/19 (k = 4.8) (a), fi = 5/13 (k = 3.6) (b), and fi = 5/7 (k = 2.4) (c). They 
are obtained by solving numerically Eqs. (51), (57), and (64), respectively, together with 
Eq. (52). The number of clusters per vertex not shown here is 1 — K + 1 and monotonically 
decreases with K from 1 to 0. 

Euler-Mascheroni constant. The sum appearing in Eq. (31) is evaluated in Ap- 
pendix A as 



= _r^--j (1 -^ + 1 -, + l r ^/ + ... ,46) 

for y(l —/j)N m ^> 1. This will be used repeatedly. We do not consider the marginal 
cases where is an integer. The sums in Eqs. (4) and (5) are evaluated using 
Ei (y) to give 

(ki) = INKPi and (k) = 2(L) /N = 2K, (47) 

in the limit N -> °° if i^P, < 1 for all j or equivalently, K < Qv(i") = d(A^ 1_ ' u ). 
Thus is {L)/N. Under the same condition, the third term (1 — e~ iPj ) in 
Eq. (6) does not contribute to (k 2 ) = (l/N)Y^\(kf) in the limit N ^ °°, which 
gives (k 2 ) = (l/N)Y!jLA(ki) + (ki) 2 ] - Moreover, rewriting Eq. (3) as 

log ft (cD) = £ log [1 - (1 -to)(l -e- 2NKFiP J)} (48) 
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and expanding the right hand side as a power series in ( 1 — to) to apply the result 
of Appendix A, we find that Eq. (44) holds for all range of < /j < 1 and K finite. 
Note that (1 - e~ 2NKPiP i) is not small when 1/2 < p < 1 but the final result is the 
same as that for < p < 1/2 where ( 1 - e - 2NKPiP i) ~ INKPiPj < 1 holds. Thus the 
degree of each vertex ki follows the Poisson distribution and the branching process 
approach and the Potts model approach are equivalent for the static model as long 
as m and s are concerned. This is because the finite clusters remain effectively trees 
for all K. 

For convenience, we divide the range of fx into the three cases, (I) < /j < 1/3, (II) 
1/3 < n < 1/2, and (III) 1/2 < /u < 1. 



5. 1 Degree distribution 



The asymptotic behavior of the degree distribution Pd(k) is related to the behav- 
ior of its generating function g(cfl) = Y,kPd(k)(8 k for (0 near 1, which is equal to 
Ei(2AT(l — (o)) in the limit N — > °°. The latter condition is necessary for the ap- 
proximation Pi(l — Pi) ~ Pi to be valid. From Eqs. (A. 2) and (A. 6), the degree 
distribution Pd(k) is given by 



c\k » for 1 < k < K 



co=0 I forfc>fc max , 



(49) 



id 



where fc max is equal to (k\) i.e., fc max = 2K(1 -fi)^, c\ = (l/fi)[2K(l -//)] 1//J and 
C2 = (1/N) Y<7=o(~kmax)%\p(k + r)]/rl ~ e~ kmdX /N. From now on, we assume that 
K e <^K<^K u with ^ = AT-^/(2(l -//)) and ^ = A^ 1 ^/(2(l -//)), for which 
1 ^max <^ N so that there exists the regime of k where the degree distribution 
follows a power law p^(fc) ~ k~^ with 



(50) 



Since we are interested in the range < fi < 1, the degree exponent X is larger than 
2. 



5.2 Giant cluster size 



The giant cluster size m can be evaluated by Eq. (33) with hi = NPiy. In terms of 
Ei (y) evaluated in Appendix A, it is simply represented as 

m=l-L l (y) (51) 
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with y obtained by solving Eq. (32) 



^- = jrP i (l-e- h <) = l+L[(y), (52) 

i=\ 



where T,[(y) = (d/dy)T.\(y). 

When y is small, Eq. (52) is expanded as 



(I) ■ y ^-rf.l 



2K 2K C 2(1 —3/i)' 



(ii) ±^-r(i-iy±^- y h\ 

y ' 2K 2K C \ uj H 



for the three ranges of (I), (II), and (III), respectively. For (I) and (II), the char- 
acteristic value K c defined by 

appears, which is just Eq. (38) with P t in Eq. (45). 

When K < K c for (I) and (II), or K = for (III), Eqs. (52) or (53) has the solution 
y = 0, and therefore, the giant cluster size is 

m = 0. (55) 

That is, there is no giant cluster for K < K c (I, II) or K = (III). But a non-zero 
solution for y occurs when K > K c (I, II) or K > (III). It leads to the following 
giant cluster size m: 



m 2(1-3//) 

(I) m ~ y ~ — — -A, 



(II) m~y 

i 

.'2(1-^ 

(III) m ~ y : 




iP^ 1 , (56) 
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for A = K/K c — 1 (I, II) or K (III) small and positive. Here the relation m ~ y comes 
from Eqs. (51) and (A. 6). The giant cluster size is finite for A finite and positive, 
and thus K c is the percolation transition point. If we define a critical exponent |3 by 
AP, its value is 1 (I) and n/{\ - 2/j) (II). For (III), m is finite for K finite, but 
K = is not a percolation transition point, which will be investigated further below. 

The giant cluster size m as a function of K, which can be obtained numerically from 
Eqs. (51) and (52), is plotted for the case of /j = 5/19 (k = 1 + = 4.8), 5/13 
(k = 3.6), and 5/7 (X = 2.4) in Fig. 2. 

We mention that in Ref. [16], some rigorous bounds of the giant cluster size are 
derived for an ensemble similar to ours, but with a different form of the probabilities 
of adding edges so that their results apply only to the case k > 3 of the static model. 



5.3 Mean cluster size 



The mean cluster size s or the susceptibility %( 1 , ho — > 0) in Eq. (34) is represented 
in terms of Ei (y) in Appendix A as 



Fi(y)] ; 
(2K)-i -m y y 



*= £ i(y) + ,^-i U,, v (57) 



where y is the solution of Eq. (52) and Tf((y) = (d/dy)L[(y). 
When K < K c for (I) and (II), y = and thus 



2KK C 

5=1 + - (58) 



K c -K 



since E(0) = -E'(0) = 1, and L'((y) = l/(2K c ). On the other hand, when K > K c , 
y is non-zero but given by Eq. (56), and one can see that E(y) ~ 1, E'(y) ~ — 1 and 



E"(v)- { 2 ^ * c (59) 
* ^_i^A (ii) 



From these relations, the mean cluster size for K around K c is obtained as 



where 



( ~ A) , (60) 

C ± (A>0) 



2K C , ( (61) 



T^^c (II)- 
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s diverges at K c both for (I) and (II). Thus, if we define % ~ |A| Y , then y = 1 for 
both (I) and (II). 

For (III), y, the solution of Eq. (52) is zero only when K is zero. We suppose that K 
is non-zero but small in the thermodynamic limit N — > °°. Then, from Eq. (A.6), it 
follows that 



J?W = ^. (62) 
where B = ( 1 - //) V/* | r( 1 - 1 ///) | . Then the mean cluster size is 



for small ^. As shown in the numerical solutions for the mean cluster size s ob- 
tained from Eqs. (57) and (52) plotted in Fig. 2, the most important feature of s 
for (III) is that it does not diverge at any value of K but has only a finite peak at 
K p 2 = 0(1). It implies that there is no phase transition for (III), i.e., 2 < X, < 3. 



5.4 Number of loops and clusters 



The number of loops per vertex (N\ oov )/N, which we denote by £, is also repre- 
sented in terms of Ei (y) as 

£ = -F (y, 0) = - 1 + y - -^y 2 + E , (y ) , (64) 
with y being the solution of Eq. (52). 

When K < K c for (I) and (II), and K = for (III), the value of y is zero and Ei (0) = 
1, which leads to 

£ = 0. (65) 

On the other hand, when K > K c for (I, II), or K > for (III), the value of £ is not 
zero. From the behavior of Ei (y) for small y and Eq. (56), one can see that for A > 
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(I, II) or K > (III), 



( 2 (1-3^ ,3 
3 (1-2^)3 A 



-TT- 



(i-2^)r(i-i 



2,-i ( 2(i-/i)^|r(i-i)| 



A"^ 



(I) , 

(II) , 

(III) . 



(66) 



The exact solutions for £ are shown in Fig. 2. The number of clusters is simply 
related to £ as (C) /N = 1 - K + £. 



6 Cluster size distribution and largest cluster size 



Beyond the largest cluster size or the mean cluster size, the whole distribution of 
cluster size P(s) for the static model can be derived from Eqs. (32) and (33), which 
gives the parametrized equations for t{z) = 1 — ra(l,^o = — lnz) as 



1 _ JL i_ JL 

1 2K 1 2K 



N 



v(z) = ±-Ye- NP » = zU(y), (67) 
iV i=\ 

where £i(y) in Appendix A is used. P(s) is obtained by P(s) = ( \ / s\)(d s / dz s )"P (z)\ z =o- 
In particular, P(s) for s 3> 1 is contributed to by such a singular term as (zo — z) x 
with x a non-integer in <P (z) . The functional form of <P (z) depends on for 1 < / < 
N. In this section, we solve Eqs. (32) and (33) when Pi is given by Eq. (45) to find 
the cluster size distribution P(s). Furthermore, we derive the largest cluster size (S) 
before a giant cluster appears through the following relation [25] 

I^) = l-^, (68) 

which is equivalent to the relation m(l,/zo — > 0) = (S)/^ in the limit ho{S) — > °° in 
Eq. (16). 



17 



6.1 n = : Erdos-Renyi model 



Before considering the case of < fi < 1 in Eq. (45), we first consider the Erdos- 
Renyi model with P l ■ = l/N corresponding to the case of fi = 0. In this case, the 
parametrized equations for <P{z), Eq. (67) is simply written as 



<B{z)=ze~ y . (69) 

If we consider y as a function of z, it has the properties y(z = 0) = 2K and dy/dz< 
for < z < zo = e 2K ^/2K. On the other hand, = 0) =0 and (p(z) is an 
increasing function of z. By substituting z = 1, we find that the giant cluster size 
m=l— fp(l)is non-zero for K > K c = 1/2, and especially, m is given by m ~ 2A 
for0< A = K/K c -l<<: 1. 

Around zo, the function y(z) becomes singular as y(z) c± 2^— l + [4,K'e 1_2A: ] 1//2 (zo — 
z) 1 / 2 . Also, (z) has the square-root singularity at zo as 

1 /e 1 - 2 ^ i 
*®~!k~\~K~) {z °- z)1 - (70) 

Differentiating 2> (z) at z = 0, one can obtain P(«), which is given for large s as 



K 



z=0 



A-2K\ 1 



r(,+i)r(-i) 



g ? 1 3 s_ 



(71) 



where r(— 1/2) = — 27C 1 / 2 is used and so = 1/lnzo. One can notice that zo = 1, 
So _> oo 5 and P(s) ~ s~ 3 / 2 at When |A| <C 1, the cut-off so is approximately 
2/A 2 . 

The presence of the cut-off so means that a cluster of size s larger than so can 
be found only with the exponentially small probability ~ e~ s / s °. Thus, the largest 
cluster is as large as sq before a giant cluster appears. (S) increase as K approaches 
K c . However, in finite size systems, the largest clusters cannot grow infinitely as 
K — > K c , which is obvious from Eq. (68). Suppose that (S) is much less than so- 
Then, one can easily see that S ~ N 2 ^ applying Eq. (71) to Eq. (68). It indicates 
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that in the regime of K where so > N 2 ^, or -1 < AN 1 / 3 < 0, (5) is o(N 2 / 3 ). For 
K > K c , the largest cluster size in the finite size system is given by Nm ~ 2NA only 
when NA 3> iV 2 / 3 . To summarize, (S) is given by 



£ (Atf*«-1), 

7V§ (|AAT3|<1), (72) 
NA (AA^ > 1). 



The regime of AT satisfying | AA^ 1 / 3 | <C 1 in finite size systems shrinks to a point 
in the thermodynamic limit N — > °°, which we will call the critical regime. If we 
introduce a scaling exponent V to describe the critical regime as lA/V 1 /^ <l,v = 3 
in the ER model. 



6.2 The case (I): < p < 1/3 



As shown in Appendix A, Ei (y) has a singular term with the /j-dependent exponent 
in its expansion in y, which allows <E (z) to have singularity other than the square- 
root one for the ER model. 

The first relation of Eq. (67) is expanded in y as, using Eq. (A. 2), 

CO 

(73) 

n=Q 

where the first few coefficients are 



a = 1, 



2K, 



1 

2K 



A 

2K' 



a 2 = - 



AKKr 



(74) 



Here, the critical point K C (N) is given by 



K C (N) = 



2(1 -^2^-1^(2^- 



(75) 



When < fi < 1/2, the solution y of Eq. (73) with z = 1 is and thus fp(l) = 1 
for a\ < while y is a positive value and rP (1) < 1 for a\ > 0. Therefore, K C (N) is 
the percolation transition point and indeed, converges to the K c in Eq. (54) in the 
thermodynamic limit N — > °°. However, when 1/2 < a/ < 1, #c(A0 = 0{N~^~ V >) 
and y satisfying Eq. (73) with z = 1 is ~ N^^A for K > K C (N), which goes to zero 
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in the thermodynamic limit. It means that K C (N) is not the percolation transition 
point for 1/2 < /a < 1. 

When y > /V - -" / ( 1 - n) , Eq. (A.6) should be used and thus 

L^— ij 

Z=1 ~J¥ + I ^+^' + -, (76) 

n=l 

where a\ = \/{2K c ), a' n = a n for n>2, and A = T(2- -/i) 1 ^" 1 . 

Similarly to the ER model, when < < 1/3, the value of m = 1 — ( 1 ) is nonzero 
for ^ > K c = 0(1), and thus K c can be identified with the percolation transition 
point. 

The leading singular term of the function y(z) varies depending on z. First, when 

1 - s m l < Z < 20 = 1 + S n 1 wim 



s = 



16^ 2 |a 2 | 

2K{\ -fi) 2 N 2 » 



A< 



4£|a 2 | 



2K\a 2 \-A(l-fi)N» V {l-n)NP 
the function y(z) is represented as, from Eq. (73), 



(77) 



y(z) 



T 1 



(78) 



\ a 2\ 



Notice that y(z) with z> 1 — s^ 1 satisfies the relation y(z) <C N ^/(l — //). Next, 
when z <C 1 — s~ l , y(z) is expanded from Eq. (76) as 



^4^b + |a2 ' 4 



zq-zA hA 



A (ZQ-Z )2 
~r 1 



r 1 



a 2 



(79) 



This implies that y(z) has the square-root singularity as in Eq. (78) except for the 
case of A < and 1 — s 1 <C z <C 1 — s^ 1 , where y(z) is given by 



/ 2K\» , i_i 



2£\ ? 



(80) 



Such regime of z exists when sq <C s„ 



Different singularities of y(z) depending on the range of z and A shown in Eqs. (78), 
(79), and (80) are inherited to <P (z) by the relation <E (z) = z£i (y), and in turn, cause 
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P(s) to behave distinctively depending on the range of s and A. When A <C — N 
sq is much less than s m and 



_3 -X 

s 2e s o (l<s< s ) 



(s>s w ). 



(81) 



P(j) for 1 < s < 5 is related to Eq. (79), that for s < 5 < s ra to Eq. (80), and that 
for 5 3> 5 m to Eq. (78). On the other hand, when A 3> —N~^, the regime of z where 
y(z) follows Eq. (80) vanishes but P(s) is simply given by 



_3 _ J_ 



(82) 



for all s. 



When A <C — 7Y the cluster size distribution decays exponentially for s 3> s m , so 
the largest cluster is as large as Sm ~ Af/|A|. Wben A 3> — N the cluster size 
distribution takes the same form as that for the ER model, and the critical regime is 
specified with the same exponent v = 3. Consequently, 



g (A « -AT*), 

(-AT"«A«-Ari), 



(|A#3| < 1), 



(5)~< 

NA (AA^i > 1), 
where Eq. (56) is used to obtain (5) for A > N~ 1/3 



(83) 



6.3 The case (II): 1/3 < p < 1/2 



Eqs. (73) and (76) are valid also for 1/3 < n < 1/2. However, one should note that 
the singular term y 1 /^ 1 is the next leading term to y in Eq. (76), which causes 
T (z) to have an /^-dependent singularity other than the square-root one even in the 
critical regime. 

Now we consider the singularity of the function y(z). As for < /j < 1/3, when 
1 - s m l < z < zo = 1 + s 1 with so and s m in Eq. (77), y(z) is given by Eq. (78). 
Note that a^ = (A 7 ^ -1 ) for fi > 1/3. y(z) for z <C 1 — s m l is obtained by inverting 
Eq. (76). When A > 0, if the regime 1 - s^ 1 < z < 1 - s m l with 
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exists, y(z) in that regime is given as 



y(z) 



2K\A\(l-fi) 



l-2fi 



+ 
+ 



2// / A \ r T 
1-2/j \2K 



(Zc-Z)- 



(85) 



with z c = 1 + s c 1 • When A < and 1 — s c 1 z 1 — s m l , y(z) is given by 



*)*£(l- l )-|A|(^J'(l-r)H + 



2£\ ? 



(86) 



Finally, both for A > and A < 0, if z is in the regime z <C 1 — max{5 c . 1 , s m l }, y(z) 
exhibits the following singularity as 



y(z) 



l-z\ l ~» 



+ ■ 



(87) 



The functional form of the cluster size distribution varies depending on s and A. 
When A < -AH 1 - 2 -"), t{z) has three different singularities, Eqs. (78), (86), and 
(87), in the corresponding regimes of z, and thus P(s) is given as 



i 



S (1<S<S C ), 

(|A|j)"^ (> c «s<s w ), (88) 

3 — -L 

s~2 e s o ( s ^>s m ). 



When —N~( l ~ 2 ^ < A< aM 1 - 2 ^ ,p( z ) i s contributed to by Eqs. (78) and (87), 
which leads to 



i 



, S (1 <5<5 m ), 

P{s)*{ _ 3 x (89) 

Finally, when A 3> N^^^ 2fJ \ the regime of z where Eq. (78) is valid disappears, 
and 

3 A ( 90 ) 

as Eqs. (85) and (87) imply. 

The largest cluster size follows s m up to A ~ — A^ 1-2 *") beyond which (S) is 
(N 1 ^) as shown by Eq. (68). The comparison of N 1 ^ and the giant cluster size 
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mN ~ NA^/t 1 2 ^ given in Eq. (56) indicates that the largest cluster size is given by 
the latter when A 3> N~( l ~ 2lJ K In summary, the largest cluster size is 



(S) 



IAI 



(AN 1 - 2 ^ -1), 
N 1 -^ (lA^V 1 -^! < 1), 
NA^ (AN 1 - 2 ^^ 1), 



(91) 



and the exponent v is l/( 1 - 2/j) for 1/3 < /u < 1/2. 



6.4 77ze case (///): 1 /2 < /j < 1 

For 1/2 < //< 1, the value of K C (N) given in Eq. (75) is o(N^^^) and thus, 
for all finite K > 0, the giant cluster size m = 1 — fp(l) is non-zero in the ther- 
modynamic limit /V — > oo. However, m is given as ~ K M ^ 2 ^~ l \ which is o(N~ >J ) 
around K C (N) and vanishes in the thermodynamic limit N — > °°. We consider finite 
size systems where K C (N) is non-zero but finite and investigate how the cluster size 
distribution and the largest cluster size behave around K C (N). 

Following the same step as for < /j < 1/2, one finds that when 1 — s m l <C z <C 
1 +s Q 1 with so and s m in Eq. (77), y(z) is given by Eq. (78). Eq. (76) applies to the 
case of z 1 — s~ l . If the regime 1 — z 1 — s m l with 



^ = 



2ju— 1 



2^-1 1-M 



(92) 



exists, y(z) is given by 



y(z) ~ 2£A 



1 — //^ ^ 
A* 



2u / 1 — t/\ 



(Zc-Z)2 + ---, 

(93) 



2//-1 V A* 

with z c = 1 + s^T 1 . In the regime z 1 — maxjs^ 1 , s m l }, y(z) is given by 

y (z) ~ 2ff( 1 - z) + A (2K) ? ( 1 - z) ? ~ 1 + • • • . (94) 



When K < K c (N) , the generating function <£ (z) is evaluated by Eq. (78) and Eq. (94) 
to give the cluster size distribution P(s) as follows: 



1 < s < s m ), 



P( s ) - f 

i 3/1-1 _3 -X 



(95) 
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Fig. 3. Plots of mN^^ vs. K for fj = 5/19 (a) and fj = 5/13 (b). The number of ver- 
tices N is 10 4 (A), 10 5 (O), 10 6 (D), and 10 7 (O)- The data cross at K c = 0.435(1) (a) 
and K c = 0.306(2) (b), respectively, which are in accordance with 0.436 and 0.305 from 
Eq. (54). 

where the factor jV _ ( 3 ^ _1 )/ 2 comes from |a2|~ 1//2 . Eq. (95) is valid when s m < s c . 
When K 3> K C (N), the regime of z where Eq. (78) applies disappears and (P(z) is 
evaluated using Eqs. (93) and (94). Consequently, P(s) is 



P(s) 



(96) 



^2(2^-1) s - le s c ( S >5 C ) : 



for K^K C (N). 



The largest cluster is as large as s m approximately up to ^(A^) and is given by mN 
in Eq. (56) for K > K C (N) . That is, 



(S) 



KN» (K<Kc(N)), 



(97) 



This result implies that the largest cluster size is still 0(N l ~ M ) even when K > 
K C (N) unless K > K C {N), consistent with Eq. (56). The conditions K < K C {N) and 
K > K C (N) can be rewritten as A < 1 and A > 1 with A = K/K C (N) - 1, and thus 
the exponent v to define the critical regime is infinity for 1/2 < /j < 1. The absence 
of a critical regime with a finite v means that there is no percolation transition in 
the evolution of scale-free graphs with 2 < X < 3, as the absence of a divergence in 
the mean cluster size does. 



7 Numerical simulations and finite size scaling 



In this section, we derive the finite-size scaling forms of the giant cluster size m, the 
mean cluster size s, and the number of loops £ and check against their numerical 
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AJV 3/13 AAf 3/13 AN 3113 

Fig. 4. Data collapse of scaled giant cluster sizes mN$' v , mean cluster sizes s/N 1 ^, and 
number of loops £N plotted vs. AN 1 ^ for fi = 5/19 (upper) and ^ = 5/13 (lower). N is 
10 4 (A), 10 5 (O), 10 6 (D), and 10 7 (Q)- 



16 I 1 1 1 j 7 




A A 



Fig. 5. Data collapse of the scaled giant cluster size (a) and number of loops (b) vs. A for 
H = 5/7. N is 10 4 (A), 10 5 (O), 10 6 (D), and 10 7 (O)- 

data from simulations for the static model. Investigating the cluster size distribution 
and the largest cluster size, we identified the critical regimes and the corresponding 
scaling variables, which enables us to predict the scaling behaviors of m, s, and t. 

Let us sketch the numerical procedure briefly. At each step, we pick two vertices 
i and j, with probabilities Pj and Pj, respectively. Then we put an edge between i 
and j unless there is one already. Repeat the procedure until there are L edges made 
in the system. The generation of random numbers with non-uniform probability 
density Pi is the most time-consuming part here. One efficient way is to use, e.g., 
Walker's algorithm [26] combined with the so-called Robin Hood method [27], 
which is explained in Appendix B. This method is exact and takes time of O (N) to 
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0.001 0.01 0.1 



K 

Fig. 6. Mean cluster size s as a function of K in semi-logarithmic scales with /u = 5/7 for 
N = 10 4 (A), 10 5 (O), 10 6 (D), and 10 7 (O)- In addition to the peak at K pl , another peak 
is shown at K p 2 ~ 0.1 iV for N = 10 5 , 10 6 , and 10 7 , respectively. The solid line represents 
the exact solution obtained from Eqs. (52) and (57). The measured values of K p \ (O) as a 
function of N are plotted in the inset together with the guide line whose slope is 1 — 2/u for 
comparison. 



set up a table and O (1) to choose a vertex, making a large-size simulation feasible. 
Once a graph is constructed, the clusters are identified by the standard breadth-first 
search, during which we can extract the size of the largest cluster, the mean cluster 
size excluding the largest one as defined in Eq. (17), and the number of clusters 
within the graph. Recall that the number of clusters is related to the number of 
loops by Eq. (18). We choose three values of ^ /u = 5/19 (k = 4.8), 5/13 (k = 3.6), 
and 5/7 (k = 2.4) as representatives of the regime (I), (II), and (III), respectively. 
We perform the simulations for the system size /V up to 10 7 , and the ensemble 
average is evaluated with at least 10 3 runs. 

Now we first consider the regime < /j < l/2orX>3. The scaling behaviors of 
the giant cluster size m in this regime are shown in Eqs. (83) and (91) and can be 
written as 

m = Arf*I' (I j I) (AA^), (98) 
where the scaling exponents (3 and v are given by 



i (o<//<i), 

3 (0 < n < \ 



2)i 



(5 <i"< 2>- 



(99) 
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The scaling function ^(iji) (x) behaves as 



„ const, (x <C 1), 
xv (x> 1). 



The critical point can be found numerically by plotting mN$l y vs. K as shown in 
Fig. 3, which are in good agreement with Eq. (54). In Fig. 4, the plots of mN$l v vs. 

are shown for ^ = 5/19 and 5/13, respectively. The data collapse confirming 
the scaling behaviors in Eq. (98). 

In the critical regime, the cluster size distribution follows the power law P(s) ~ 

s~ x with x = 3/2 for (I) and 1/(1 — /j) for (II). Then the mean cluster size in the 

critical regime is given by s ~ (S) 2 ~ x since s = L^{s) sP(s). Since (S) ~ A^ 1_ P/ V 

in the critical regime as in Eq. (98), the mean cluster size is represented as s ~ 
^(i-p/v)(2-x) = N i/v^ which leads to> combined with Eq ( 60 ), 

s = N** m {AN*), (101) 
where the scaling function ^(i^i) (x) behaves as 

„ const. (x<£i 1), 
*fcn)to~< , (102) 
x- 1 (x>l). 

Plots of s/N 1 /* vs. AN 1 /* for /j = 5/19 and 5/13 are shown in Fig. 4. 

Since the giant cluster size and the mean cluster size correspond to the first and 
second derivative of the free energy with respect to the external field, respectively, it 
is natural that both have the same scaling variable. Therefore, the number of loops, 
corresponding to the free energy itself, is also described in terms of the scaling 
variable AiV 1 ^ with v in Eq. (99). Considering the number of loops for K > K c 
given in Eq. (66), one can see that 

£=^a m (AN^), (103) 
where the scaling function Q( L n)(X) behaves as 

, const. (x< f 

. (104) 

y (x»i). 

This implies that £ is o(N~ l ) in the critical regime, which is supported by the 
numerical simulation results in Fig. 4. 
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Next, we consider the regime 1/2 < fi < 1. From the largest cluster size shown in 
Eq. (97), the giant cluster size m can be written as 



m~N-» ¥ (m) (A), 



(105) 



where the function ^(m) (x) behaves as 




const, (x <C 1), 
x^~ l (x>l). 



(106) 



Notice that m increases smoothly passing K C (N) as manifested by A not scaling with 
N. Similarly, the number of loops i represented in terms of the scaling variable A 
around K C (N) should exhibit the following scaling behaviors to satisfy Eq. (66): 



Data collapse of mN M and £N vs. A is shown in Fig. 5. 

The numerical data of the mean cluster size are shown in Fig. 6. As /V increases, 
the mean cluster size s approaches the exact solution in Eq. (57) represented by the 
solid line in Fig. 6. It does not diverge at any value of K, but instead its peak height 
decreases as N increases. We additionally find that the mean cluster size s has a 
small peak at K p \, which scales as N l ~ 2fJ as shown in the inset of Fig. 6. The value 
of K p i is not equal to K C (N) although they are the same order of N. The reason for 
the peak at K p \ is as follows. As K increases, the largest cluster size (S) and the 
mean cluster size s = L*^(5) sP(s) also increase. However, as K approaches K C (N), 
the cluster size distribution P(s) begins to develop the exponential decaying part in 
its tail, i.e., for s 3> s c . s c decreases with increasing K after K passes K c . At K p \, 
(S) and s c are equal. After K passing K p \, the mean cluster size is dominated by 
s c , which makes s decrease for K > K p \. However, the mean cluster size increases 
again as soon as K becomes much larger than K p \ or K C (N) because the prefactor 
K 1 ^ of the cluster size distribution P(s) for 1 s <S s c increases with increasing K. 
The mean cluster size decreases only after the second peak at K p 2 = 0(1), where 
s c = 0(1), as shown in Fig. 6 as well as in the exact solution in Fig. 2. Since 
K p 2 = (1), a giant cluster exists around K p 2. 



(107) 



where the scaling function £2( ni )(x) behaves as 




(108) 
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8 Summary and Discussion 



In this paper, we have studied the percolation transition of the SF random graphs 
constructed by attaching edges with probability proportional to the products of two 
vertex weights. By utilizing the Potts model representation, the giant cluster size, 
the mean cluster size, and the numbers of loops and clusters are obtained from the 
Potts model free energy in the thermodynamic limit. Our general formula for the gi- 
ant cluster size and the mean cluster size are equivalent to those results obtained for 
a given degree sequence if the latter expressions are averaged over the grandcanon- 
ical ensemble. The Potts model formulation allows one to derive other quantities 
such as the number of loops easily. Using this approach, we then investigated the 
critical behaviors of the SF network realized by the static model in detail. Further- 
more, to derive the finite size scaling properties of the phase transition, the cluster 
size distribution and the largest cluster size in finite size systems are also obtained 
and used. We found that there is a percolation transition for X = 1 + > 3 so that 
a giant cluster appears abruptly when K = (L)/N is equal to K c given by Eq. (54) 
while such a giant cluster is generated gradually without a transition for 2 < X < 3. 
Thus the process of formation of the giant cluster for the case of 2 < X < 3 is 
fundamentally different from that of X > 3. 

In particular, the difference between the SF graphs with X > 3 and 2 < X < 3 is man- 
ifested in the mean cluster size. For X > 3, as K or the number of edges increases, 
many small clusters grow by attaching edges, which continues up to K = K c , and 
then a giant cluster forms by the abrupt coalescence of those small clusters as shown 
in Fig. 7. Since we do not count the giant cluster in calculating (s), s decreases 
rapidly as K passes K c . Thus the mean cluster size exhibits a peak at K = K c , which 
diverges in the thermodynamic limit N — > °°. On the other hand, for 2 < X < 3, the 
role of K c is replaced by K p \ which is ~ o(N~^~^^^'~ l ' > ), but after a small peak, 
the mean cluster size increases again after passing K p \ as seen in Fig. 6: Edges 
newly introduced either create new clusters of size larger than 1 or merge small 
clusters to the larger one with its size not as large as 0(N). Only when K reaches 
K p 2, the network is dense enough for the giant cluster to swallow up other clusters 
to reduce s on adding more edges. The evolution of the static model as the density 
of edges increases is summarized in Fig. 8. 

Due to the characteristics of the power-law degree distribution of scale-free graphs, 
the cluster size distributions exhibit crossover behaviors such that they follow dif- 
ferent power laws depending on the cluster size for given X. In particular, for 
3 < X < 4, the scaling exponents |3 and v depend on the degree exponent X con- 
tinuously while they are equal to the conventional mean-field values for X > 4. We 
have not considered the marginal cases where thus X is an integer value where 
logarithmic corrections appear as seen in Eq. (A.7). 

The static model considered in this paper, though algebraically tractable, does not 
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Fig. 7. Schematic picture for the comparison of cluster evolution between the case of X > 3 
(a) and of 2 < l< 3 (b). 
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Fig. 8. Schematic phase diagram of the static model. The SF behavior of the degree distri- 
bution appears between K( and K u . A giant cluster emerges at K c for (1,11), and around K p 2 
for (III). The quantities in [• • • ] are only for (III). 

have correlations between edges that are important in real world networks. Exten- 
sion of this work to the cases allowing correlations is left for the future work. 
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Appendix A The evaluation of Ei (y) 



We derive the sums 



(A.1) 



i=l 



where P t = i^/^in) (0 < ji< 1) with = Z,f =i r» as in Eq. (45). 

Let us introduce y^ = ^^y.y> which is equal to (1 — jS)y in the thermodynamic limit 
]V^«we are interested in. Then Ei (y) = £jLi s (0 wim s(x) = g-wW^/M 
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First consider the regime where y^N** <C 1. Then Ei (y) is simply expanded as 



iii 

with [xj being the greatest integer less than or equal to x. 

On the other hand, in the limit y^/N M — > °°, the second summation gives rise to 
a non-analytic term. Since s(x) and its derivatives at x = 1 and x = N have the 
properties = otA^V^") and jW(iV) = o(Ar B_1 ) with n > 0, the 

Euler-Maclaurin formula [28] enables us to evaluate Ei (y) in the limit N — > °° and 
ywiV 7 ' — > 00 as follows: 



Ei(y)= / dxs{x)= dx-e- yN ^> 
Ji Ji N 



1 1 f 00 
= -y£ / <fc "e z 
" J yN 

1 W 1 

= -W — ,yw 
where the incomplete Gamma function r(s,.x:) is defined as 

T(s,x)= / dtf^e-t. 

Jx 

The incomplete Gamma function can be expressed as 



(A.3) 



(A.4) 



\s,x)= f 

J X 

Jo 

fX 

- \ dtf- 1 

Jo 

J x 



L--vJ 



H) n + L ^ J H) n 



w=0 «=0 



(-0" 



w=L-jJ+i 



(-ty 



n=|_-ij+l 



(A.5) 
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where it is used that T(s) = /J°rff/ s-1 [e - ' - £ n Jo (-t) n /n\], and therefore, it fol- 
lows that 



f — 11" / 



y N - 



(A.6) 



As jU — > 1/m with m an integer, a logarithmic term is developed as 



+ 



Tl n\{\ —un)~ 

(-1)' 



(m-1) 



hrm-YM + l + ^H ^-r 

2 m— 1 



(A.7) 



since 



lim 



(-i) OT „ t l n i 

m!(l— jjm)^ N \ n)^ N 



(-ysT 



(1 — jjm)T(m+ 1) 



1- 



(m-1)! 



l + (±-ra) lnvyv 
1 

_i 

2 



r(m) ^/i 



m — \ 



where it is used that r(z)T(l — z) =7t/ sinrczand r'(m)/r(m) 
\/{n + m-\)}. 



(A.8) 



-YM+L,T =1 [i/n- 



Appendix B Walker algorithm 



Suppose we want to choose discrete random numbers %i,X2, . . . ,jcat, with probabili- 
ties p\ ,p2, ■ ■ ■ , Pat, respectively, where p ; 's are arbitrary yet properly normalized as 
Yli=\Pi = 1- Walker algorithm enables one to chooser's with appropriate frequen- 
cies by picking a single random number < r < 1. To do this, however, we have to 
set up a table {(qi,yi)}. The table is made in the following way. 

(1) Initialize qfs as qi — Npi, (i = 1,2,.. . , N). 

(2) Divide xfs into the poor (qi < 1) and the rich (qt > 1). 

(3) Pick a poor, say p, and a rich, say r. 

(4) Fill the shortage of p, 1 — q p , from r. 

(5) p records the donator r as y p = r. 
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(6) q t is updated as q r <— q r — (1 — q p ). 

(7) If r becomes poor by the donation, i.e.,q r < 1, r enters the list of the poor out 
of that of the rich. 

(8) Repeat (3)-(7) until there are no poor left (hence the name Robin Hood 
method). 



The donate-and-fiH steps (3)-(7) are performed at most N — 1 times, thus the table- 
making procedure takes time of O (N). In its original introduction [26], Walker pro- 
posed the step (3) as "pick the poorest and the richest," which involves additional 
sorting operation, increasing the computational cost. The present scheme follows 
the implementation of Zaman [27]. The table-making process can be visualized for 
a simple N = 3 case as in Fig. B.l. 



Pi 



(a) 



P 3 



(b) 



Fig. B.l. Illustration of the table-making procedure in the Walker algorithm for N = 3. 
heights of the white bars in (c) indicates qi/N and the figures in the shaded boxes y,\ 



The 



With the table {{qi,yi)} at hand, one picks a random number < r < 1. Divide 
x = rN + 1 into the integer part n = [x\ and the remainder d = x — n. If d < q n , we 
choose x n , otherwise choose y n . With this scheme one can draw discrete random 
numbers xfs with appropriate probabilities p/s. 
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